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INTRODUCTION 


An important goal of geodesy is to determine the anomalous potential 
and its derivatives outside of the earth. Representing the surface 
anomalies by a series of spherical harmonics is useful since it is then 
possible to do a term by term solution of Laplace’s equation and upward 
continuation. This paper addresses the problem of finding such a 
spherical harmonic series for anomaly values given on an equiangular 
surface grid. (This is a first step toward the more complicated problem 
of finding a function such that locally averaged values fit a grid of mean 
anomalies.) Three approaches to this fitting problem are discussed and 
compared: the discrete Fourier technique, the discrete integral 

technique, and a new approach by this author. The peculiar nature of the 
equiangular grid, with its increasing density of (noisy) data toward the 
poles, causes each method to exhibit a different type of difficulty. The 
new method is shown to be practical as well as precise since the numerical 
conditioning problems which appear can be successfully handled by such 
well-known techniques as a (simple) Kalman filter. 

DISCRETE FOURIER METHOD 


The discrete Fourier method [ Dilts , 1985] uses a discrete Fourier 
series to represent both the longitude and latitude variation of the 
desired function. The data at the (i, j ) grid point on a grid of N 
latitude and 2N longitude intervals can be uniquely represented by the 
double Fourier series. 
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The discrete Fourier method makes its modeling assumption at this point by 
choosing the function off the grid points to be given by this same double 
Fourier expansion. Comparison of the continuous spherical representation 
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and expansion of the normalized Legendre polynomials 
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to the function modeled as in 
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Equation (1) then yields 

A for |q|=N 
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0 otherwise. 
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A solution exists for the upper limit L, equal to infinity. It can be 
expressed as 
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where the "inverse" coefficients are obtained from 
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for 0 between zero and 2 71 radians. 

The shortcoming of this approach is the need for an infinite number 
of terms to solve Equation (4) for arbitrary A qm (representing the data). 

Small amounts of noise in A qm can lead to the presence of terms in the 

double Fourier expansion (Eq. (D) which are not present in the gravity 
field and which have infinite derivatives at the poles. Truncation of the 
series is the strategy for coping with this difficulty After truncation, 
the function will no longer match the gndded data, and the degree of 
discrepancy is not under the analyst s control. 

DISCRETE INTEGRAL METHOD 

The discrete integral approach has been widely used (see for example 
Colombo [1981]). It approximates the continuous inversion integra or 
the spherical coefficients by a discrete, weighted sum. 
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The weights W. are usually chosen to be the grid block areas. The 
difficulty with this approach is that the discrete P nm <V are not 

orthogonal on the equiangular grid. As a result, aliasing occurs, 
and the resultant spherical expansion does not match the grx e a * 

The expansion is truncated at degree N or less, and the amount of the 
discrepancy is thus onlv indirectly under the analyst s control. 

Comparison with the preceding technique is obtained by using the expansion 
of Equation (6) in the above expression (with the weights proportional 
area and the interval extended to 2^): 
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Comparison with Equation (5) shows that it corresponds to the leading term 
in the above expression. Thus, taking the degree of the discrete integral 
expansion to infinity does not appear to reproduce the gridded data. 


NEW METHOD 


The third method is newly presented by this author. It uses the 
Fourier representation of the data (Equation (1)) but makes its modeling 
assumption in the spherical domain. Comparison to the spherical expansion 
only at the grid points yields 
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This differs from Equation (4) since it is the result of a discrete 
comparison at the grid points (using the periodic nature of the discrete 
exponential) and not a comparison of continuous functions. If L is chosen 
equal to N+|m|-2 (except L=N for m=0) , Equation (9) then becomes an 
invertible matrix equation (with E indicating the sum of the B terms): 

EC = A and then C = E _1 A. (10) 

Since the inverse yields a precise fit at the data points, the modeling 
assumption is that all the C T s are zero for n greater than N+|m|-2. The 
continuous function resulting from using these C’s in a spherical 
expansion thus reproduces the data and has a finite number of terms. 

Since L<2N the elements of the matrix E are easy to compute: at most two 

of the B terms in Equation (9) are non-zero. Even for terms of degree 
less than N, this solution is different from the discrete Fourier case, 
Equation (5), since (ZE) is not the identity. 

The difficulty with this method is that the matrix E becomes 
ill-conditioned for large values of the order m. There are, however, many 
well-known and trustworthy techniques for dealing with such problems. A 
few such techniques are summarized below. 

o Perform the transformation of E to the identity in a column by 
column fashion, stopping when the conditioning becomes a 
problem. If this process is stopped at the column for degree N, 
the discrete Fourier approximation is obtained. Further steps 

toward finding E ^ constitute improved approximations. 
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Invert the matrix <E+6l) for a small 6 and use it instead of 


o Use a simple Kalman filter 

A = EC + V; C = E T (EE T + yI) _1 A (11) 

where the measurement noise, V, has variance y^I and the prior 
uncertainty on C is y 2 I. Then y = and is a small 

quantity. 

o Use a more complicated Kalman filter with detailed models for 
the noise and for the initial uncertainty. 

All of these strategies yield results which are not overly sensitive to 
noise. By adjusting the parameters in these methods, the analyst can 
control how close the reconstructed function comes to the gridded data 
(allowing only for small deviations consistent with the noise model). Use 
of the Kalman filters also has the advantage of providing uncertainties in 
the estimated spherical coef f icients . 

SUMMARY 


The problem of fitting a smooth function to data given on an 
equiangular spherical grid has been discussed. Two existing approaches 
were summarized and a new approach was presented. Each approac wa ^ oun 
to possess an area of difficulty resulting from the properties of the 
equiangular grid. Well-known techniques (such as Kalman filtering) are 
available as practical strategies for dealing with the numerical 
conditioning in the new method. As a result, the new method is practical 
and capable of reproducing the gridded data to a precision consistent wit 

the noise model. 
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